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ABSTRACT 

We study explosions of stellar models using a one-dimensional Lagrangian hydrody- 
namics code. We calculate how much mass is liberated as a function of the energy 
of explosion for a variety of pre-explosion polytropic structures and for equations of 
state with a range of radiation-to-gas pressure ratios. The results show that simple 
assumptions about the amount of mass lost in an explosion can be quite inaccurate, 
and that even one-dimensional stellar models exhibit a rich phenomenology The mass 
loss fraction rises from about 50 to 100 per cent as a function of the explosion energy 
in an approximately discontinuous manner. Combining our results with those of other, 
more realistic models, we suggest that Nova Scorpii (J1655-40) may have experienced 
significant mass fallback because the explosion energy was less than the critical value. 
We infer that the original progenitor was less than twice the mass of today's remnant. 

Key words: supernovae: general, computational; hydrodynamics 



1 INTRODUCTION 



A fundamental question in the study of supernovae is the 
fate of a star subject to an explosion of a given strength: 
is the star completely disrupted, and, if not, how much of 
the star is lost and what is the configuration of the matter 
that remains bound? Many researchers have addressed this 
question for specific cases of interest using detailed numer- 
ical simulations. To our knowledge, a precise quantitative 
relationship between the strength of the explosion and the 
fate of the outer layers has not been given before, even for 
highly idealized stellar models. The pote ntial utility of such 
a relat ionship is evident in the analysis of lFrver fc Kaloeeral 
teOOll) . where a simple "rule of thumb," introduced to esti- 
mate the amount of mass left bound in a supernova, allows 
a determination of which high mass stars leave behind neu- 
tron stars and which ones yield black holes. The "rule of 
thumb" stipulates that a portion (between 30 and 50 per 
cent) of the explosion energy is effective in directly unbind- 
ing the outermost layer s of a star; this es t imate is based on 
detailed simulations bv lMacFadven et alj i20Q]J) for a set of 
specific stellar progenitors. 

Our goal is to improve the understanding of the disrup- 
tion process by carrying out hydrodynamical calculations of 
simple, polytropic stellar models with a range of explosion 
strengths, polytropic indices, and equations of state. We in- 
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tend this sort of calculation to complement, not replace, the 
realistic, detailed simulations that are the current state-of- 
the-art in this field. Physical complications such as density 
jumps, neutrino transport, and aspherical motions are ab- 
sent from our calculations. Instead, our goal is to incorpo- 
rate the essential physics - hydrodynamics and gravity - in 
models that are easy to compute and useful to the study of 
supernovae in the same way that the polytrope itself is use- 
ful to stellar modeling. We note that this subject has been 
treated before, in a very different per turbative calculation 
iNadezhin fe Frank-Kamenetskiill963l) . Already, our results 
yield improved versions of the "rule of thumb," which we 
provide in a simple, easily applied, empirical form. Although 
a host of significant core-collapse modeling uncertainties re- 
main (hydrodynamic motions in the core, distribution of 
angular momentum within the collapsing object, neutrino- 
matter coupling, etc.), our simplified treatment represents 
an improvement in the determination of the fate of central 
remnants - providing a convenient bridge between estimates 
of mass loss based upon simple physical assumptions and 
more sophisticated, realistic simulations. 

There is considerable evidence that a supernova explo- 
sion occurred in J1655-40: the atmosphere of its compan- 
ion is contaminated with elements thou ght to be formed 
only in supernovae (Israelia n et alJll999T) . and it is likely 
that the black hole progenitor was c onsiderably more mas- 
sive than the remna nt we see today llOrosz fc Bailvnlll997l : 
IShahbaz et al.l fl999''l. There is also some evidence that the 
J1655-40 system could have remained bound only if it re- 
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ceived a substantial kick during: or short ly after the forma- 
tion of its black hole jMirabel et all2 002 ) . While our current 
models are too simple to provide definitive results for any 
particular system, we believe that the methods employed 
here suggest that the progenitor mass was less than twice 
the mass of today's remnant. Mass fallback may trigger the 
collapse to a black hole as well as pollute the companion's 
atmosphere. 

In section 2, we describe the physical set up, while sec- 
tion 3 describes the numerical code. In section 4, we give 
more detailed results and discuss how the numerical data 
were analyzed, and in section 5 we comment on our results' 
lack of dependence on the inner boundary condition. 



2 PROBLEM AND PARAMETER RANGES 

We model the supernova as a spherically symmetric explo- 
sion in a stellar model that is initially in hydrostatic equilib- 
rium. The pre-explosion stellar structure is polytropic. We 
deposit the full energy of the explosion in a small region near 
the centre of the polytrope. Using a finite-difference code we 
calculate the hydrodynamical evolution. A shock propagates 
towards the surface and the outer layers of the stellar model 
may be ejected. If it is not completely destroyed, part of 
the model remains gravitationally bound and we follow the 
evolution long enough to make an accurate estimate of the 
mass of the remnant. 

We considered a range of initial stellar structures. We 
varied the polytropic index V where P oc /. The Lane- 
Emden equation prescribes the run of density and pressure 
in the initial model; our choices for n = 1/(T — 1) span 3/2 < 
n < 4. As is well-known, the polytrope's ratio of central to 
mean density increases as n varies from to 5. This range 
subsumes typical main-sequence profiles and extended red- 
giant structures. 

We considered two equation of state treatments: ideal 
gas pressure ("EOS M": P = Pmatter with a fixed ratio of 
specific heats 7) and a mixture of gas plus radiation in ther- 
mal equilibrium ("EOS MR": P = P rad + Pmatter)- EOS 
M is suitable for stellar models of low mass (dominated 
by particle pressure at their centres) and weak explosions 
(such that the post-shock gas is not radiation dominated); 
EOS MR is needed if there is significant radiation pressure. 
We infer the temperature profile from the appropriate EOS 
and the Lane-Emden pressure-density profile. For EOS MR, 
we chose to limit ourselves to convectively stable polytropic 
models. We will discuss the condition for stability in §4.2. To 
facilitate the description of our problem's parameter space, 
let s c = P ra d(r = 0)/ Pmatter (r = 0) . As we will show in §4.1, 
for a given polytropic index n, the choice of s c uniquely fixes 
the mass of the resulting stellar model. We define our mass 
scale to be M sca ie = MQ(m p /2fi) 2 . We then chose the di- 
mensionless masses of the stellar models we exploded to be 
m = M a tar/M 3ca i e = 10, 100, 1000 for a range of (n, s c ) 
pairs, taking care to remain in the convectively stable region 
of parameter space. See TabQand FigQ 

We considered a variety of explosion energies. Given 
our interest in studying explosions which only partially un- 
bind the stellar models, we typically considered blasts with 
0.1 < -Ebiast /-Bbind < 1-5, i.e. energies of the same order of 
magnitude as a simple dimensional estimate for unbinding. 



In brief, the results we obtained were as follows. 

• The amount of mass lost as a function of explosion en- 
ergy makes a discrete jump from approximately 50 to 100 
per cent in all models, suggesting a point of instability. 

• Explosions that do not totally disrupt a stellar model 
give rise to mass loss curves that, in most cases, scale 
quadratically with explosion energy. 

These simulations did not include any sort of compact 
object at the core of our explosions. We explored the sensi- 
tivity of the mass loss to our treatment of the inner boundary 
condition. We found that to a large extent the results are 
unchanged: 

• When two extreme inner boundary conditions - a 'vac- 
uum cleaner' core that sucks up any incident material, and 
its opposite, a hard, reflecting shell - were compared, the 
mass loss results were virtually identical. 

Readers primarily interested in further details regarding 
these results are encouraged to skip to §4. 



3 THE CODE AND NUMERICAL TESTS 

3.1 Equations 

We use the inviscid fluid equations which describe mass, 
momentum and energy conservation. All calculations are 
one-dimensional with either a plane-parallel (for testing) 
or spherical (for testing and simulations) geometry. We ad- 
vance the fluid state using a finite difference approximation 
to the fluid equations (Lax-Wendroff, explicitly di fferenced, 
1-D Lagrangian code IRichtmver fc Morton! Il967| '). Shocks 
are handled with the addition of artificial viscosity. We solve 
Poisson's equation to determine the gravitational forces at 
each time step. Details are provided in Appendices A - C. 

3.2 Tests of hydrodynamics 

We tested the purely hydrodynamic capabilities of the code 
(no gravity) via comparison with the Sod shock tube (plane- 
parallel geometry) and Sedov blast (spherical geometry) so- 
lutions. For the Sod test with 7 = 7/5 (as well as for a 
range of other 7's), EOS M, various overpressures (P2/P1 = 
10, 100, 1000), and 1200 zones, we found essentially perfect 
agreement between the numerical and analytic solutions, ex- 
cept for the shock smearing over ~ 5 — 8 zones. 

For the Sedov probl em, we re-derived the solu tion given 
in Fluid Mechanics bv lLandau fe Lifshitzl Jl987ft . thereby 
finding the correction to that solution's typographical error 
(in an exponent) mentioned in Shu's Gas Dynamics. The 
correction is recorded in Appendix D. We carried out a num- 
ber of blast wave simulations, varying our choices of EOS 
and 7. For flows dominated by particle pressure we com- 
pared numerical solutions (7 = 5/3 and 7/5 for EOS M) 
with the analytic similarity solution; for flows dominated by 
radiation pressure we compared several different radiation- 
dominated numerical solutions (7 = 5/3, EOS MR) to the 
7 = 4/3 similarity solution. The radiation-dominated nu- 
merical solutions were generated by explosions yielding high- 
Mach-number shocks. A range of initial radiation-to-matter 
pressure ratios s and explosion energies were considered. 
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One simulation with large constant s ~ 1000 and relatively 
small explosion energy and another simulation with small 
constant s ~ 0.1 and large energy both yielded a radiation- 
dominated post-shock flows. 

In all calculations, the explosion was allowed to expand 
to well over 100 times the size of the initial "bomb zone." 
Comparisons of EOS M runs with analytic solutions were 
possible throughout the simulation; comparisons of EOS MR 
runs with the analytic radiation-dominated 7 = 4/3 similar- 
ity solution were meaningful only for the part of the simula- 
tion in which radiation pressure dominated matter pressure, 
approximately 4-5 expansion times. With 800 zones, the Se- 
dov test gave close (2 — 3 per cent) agreement in the relative 
density, velocity and pressure of the numerical solution and 
the analytic similarity solution for both particle pressure 
dominated and radiation pressure dominated flows except 
in the central-most region. 

Two factors contribute to the discrepancies at the cen- 
tre. First, the innermost zone was treated as an adiabatic 
expanding/contracting bubble. The entropy of this zone was 
incorrect but its mass was so small that its impact on the 
rest of the solution was inconsequential. The explosion re- 
sults were found to be almost entirely insensitive to alter- 
native methods of treating this innermost zone, provided 
the treatments were energy conserving. Second, explosive 
energy was injected in a small, but non-negligible central 
region (typically the inner 5 per cent of the mass). Quan- 
titative differences between the analytic similarity solution 
and the numerical solution occurred in the part of the grid 
used as the "bomb zone" and persisted through the simu- 
lation. These differences were not unexpected, as the point- 
like nature of the explosion in the similarity solution cannot 
be realized in any finite simulation. We also compared two 
models for the energy injection at the centre. In one, the 
"thermal bomb," an excess of thermal energy equal to the 
desired explosion energy was added by hand to the core (in- 
ner 5 per cent of the mass) of the stellar model, essentially 
creating an out-of-equilibrium hot core that then expanded 
rapidly into the stellar envelope. In the other, the "kinetic 
bomb," a linear velocity profile carrying the same amount 
of energy was added to the inner 5 per cent of the mass. 
Both methods produced identical results outside the "bomb 
zone." 



3.3 Tests of hydrostatics 

With the inclusion of self- gravity forces, we verified that 
Runge-Kutta integration of the Lane-Emden equations 
yielded stationary, stable configurations for our time- 
dependent hydrodynamic evolution equations (finite dif- 
ference scheme). We checked the long-lived stability for 
all polytropic indices and radiation-to-gas pressure ratios 
adopted in this study. Likewise, we tested that the virial 
theorem was satisfied by the initial configurations. 

3.4 Tests of self-gravitating explosions 

In the actual runs of the problem of interest, we further 
verified that the treatment of the central zone made no dis- 
cernible difference, that variations in the size of the "bomb 
zone" (3-10 per cent, for instance), caused only very slight 



(< 5 per cent) changes to the amount of mass lost in the 
explosions. We also verified that energy conservation was 
satisfied (to < 5 per cent). 



4 RESULTS 

We adopted polytropes for the initial stellar structure with 
P = kp F . The density and pressure profiles were determined 
by solving the Lane-Emden equation with the total mass 
and radius scaled to unity. We refer to this as the dimen- 
sionless solution; it depends only upon F. The dimensionless 
density-pressure distributions are the forms used in our com- 
putations. All results are likewise reported using dimension- 
less quantities (explosion energy in terms of binding energy, 
mass loss in terms of the total mass, etc.). 

4.1 Scaling of polytropes 

Let us first review the scaling of the initial polytropic solu- 
tion. For given k and r = l + l/nin the pressure-density 
relation, it is possible to generate a one-parameter family of 
scaled solutions with 

M 2-r R 3T-4 = ( fc / G )j(r) ) (!) 

with /(r) a dimensionless number depending on polytropic 
index. We can construct a polytropic progenitor without ra- 
diation pressure, but that is an idealization that is approx- 
imately correct only in the limit of a low ratio of radiation 
to gas pressure. If the matter has an ideal gas equation of 
state, P = pkT/p, there must be nonzero temperature in- 
side the stellar model, and hence nonzero radiation pressure 

Praxi — 

aT 4 . Under some circumstances, P ra d will be low 
both in the progenitor and in the ejecta after the stellar 
model explodes. The explosions of such models have a uni- 
versal mass loss fraction as a function of F and explosion 
energy in units of the stellar binding energy. 

Imposing a fixed value of the radiation-to-gas pressure 
s(r — 0) = s c at the centre of the stellar model before the 
explosion reduces the scaling of the dimensionless solution. 
These models also have a universal mass loss fraction as 
a function of F, s c , and explosion energy. The reason the 
scaling is reduced is because specifying s c determines the 
stellar mass for a given value of F independent of k: 

M = m(F)M C h4 /2 (1 + sc) 3/2 , (2) 

where M C h = (?ic/G*) 3/ V~ 2 = lM(m p / pf M Q = 
7.44M scaie , and m(T) = (45/7r 2 ) 1 / 2 [/ p (r)] 2 [/ p (r)]- 3 / 2 , with 
/p(r) = P C R A /GM 2 and f p (T) = p c R 3 /M, which are 
both dimensionless functions of polytropic index only. For 
n = 3/2, f p = 1.430 and f p = 0.7702, so m(5/3) = 6.460; 
for n = 3, f p = 12.94, f p = 11.05, so m(4/3) = 9.734. In 
the low mass limit, this implies s c tx M 2 . Thus, from Eq. 
0, the combination kR 4 ~ 3T is determined given F and s c . 
Recovering the s c — > limit is subtle, since it also implies 
low mass M. To summarize: for EOS MR, we use Eq. @ to 
solve for s c given M in terms of M aca i e = (m p /2p) 2 Mq. 

In this paper we will adopt the point of view that k is 
not known a priori and we will allow scaling of the polytropic 
solution to arbitrary M and R in cases with no radiation 
pressure. In cases with radiation pressure, although M is 
determined from Eq. , some scaling remains since Eq. 
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relates R and fc, given M and T, but does not determine 
either one separately. 

4.2 Convective Stability of Polytropes 

If we begin with the First Law of Thermodynamics, 

TdS = dE + PdV, (3) 

and use EOS MR, together with some of the relations de- 
rived in Appendix B, we can put it in the form 



TdS = f(s,n)dP/p. 



(4) 



Since dP/dr < 0, local convective stability requires dS/dP < 
0. The condition is 



n > rtcrit = 



3(8s r + l)(s r + 1) 



(•>) 



where we have used s r — s(r) = P ra d{r) / ' Pmatter{r) to dis- 
tinguish between this function of r and the constant param- 
eter s c - For stability, we require that the local condition be 
satisfied throughout the model. The variation of s r depends 
upon n: 

p3 3 _ n 

S r (l + S r ) OC — OC P !+" . 

p 

There are three characteristic cases 



(6) 



(i) For n < 3, the largest value of s r is s c = s(r — 0), 
so the largest value of n cr it is n cr it(s c ) and therefore if n > 
n C rit(s c ), the stellar model is stable. 

• For all < s r < oo, Eq. (5) implies n cr i t (s r ) > 3/2, 
so no polytrope with n < 3/2 can be stable. 

• For 3/2 < n < 3, there is a maximum value of s c for 
which the stellar model is stable, determined from n — 
n cr it(s c ), which rises from s c = at n = 3/2 to s c — > oo 
for n — 3. For instance, at n — 2.0, s^ IAX = 0.3; for 
n = 2.5, sf Ax = 1.674 (see Tab.[T]and Fig.[TJ. 

(ii) For n = 3, s r is constant. Since n C rit(s r ) < 3 for any 
finite s r , all models are stable. 

(iii) For n > 3, s, - — » oo at the polytrope's edge (where 
P — * 0), so for any choice of s c the largest value of 
ncrit(sr) = 3. All models are stable. 

It is important to emphasize that the convective-stability 
condition described here is not an absolute stability crite- 
rion. Stellar models outside of this mass range can certainly 
exist; they will simply have convection zones. The unstable 
region in Fig. arises because of our insistence on an ex- 
act polytropic pressure-density relation (P — np T ) and our 
use of EOS MR. We restrict our studies to stable models. 
These already span a wide range of density profiles - which 
is the basic physical property that governs shock propaga- 
tion - from the highly centrally-condensed n = 4 polytropes 
to the diffuse n — 2 models. An unstable (polytropic) model 
would represent an unphysical and unnatural initial state. 
Actual convective stars involve physics beyond the scope of 
our simple, idealized, one-dimensional approach. 



4.3 Description of Analysis 

The chief way in which we shall summarize the results of an 
explosion is in terms of the mass ejected as a function of ex- 
plosion energy. We begin by discussing how we extracted the 
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Figure 1. A plot summarizing the results of the stability analy- 
sis described in §4.2. Convectively stable and unstable regions of 
the n - s c parameter space are labeled, and the models we have 
selected from this parameter space are marked. 



ejected mass from the numerical simulations. For EOS M (no 
radiation pressure) , the code was run until the remnant core 
had become stationary and had nearly reassumed hydro- 
static equilibrium, i.e., it had local gas velocities near zero 
(< 10~ 7 i? s teHar/dynamical time) and satisfied the virial the- 
orem. The mass loss was determined by finding the location 
in the Lagrangian grid of the outermost outgoing zone for 
which the local energy per unit mass (sum of kinetic, ther- 
mal and gravitational contributions) changed from negative 
to positive (see arrow in Fig. . A graph of the local en- 
ergy density for a typical stellar model after an explosion is 
included as Fig. [5] 

A drawback of this method is apparent in Fig. [5] 
Though there are distinct portions of the stellar model that 
can definitely be said to be either remnant or ejecta, there 
is also a small region with nearly zero energy, resembling 
an atmosphere around the remnant. These atmospheres did 
not appear in all explosions - typically, they occurred when 
-Ebiast ~ 0.7 - 1.0 -Ebind- In some cases, it was adequate sim- 




Figure 2. The asymptotic dimensionless local energy per unit 
mass (U / (G M s tar / Rstar)) (solid line) for n = 3 polytrope (no 
radiation pressure) and iSblast = -^bind after the central core has 
reattained hydrostatic equilibrium. The dash-dotted line gives the 
dimensionless local velocity (»/(GM,t ar /Rstar) 1 ^ 2 )- The dashed 
line is a reference line for zero energy and velocity. 



ply to run these models longer, with a clear bifurcation point 
eventually emerging. 

When we began to us EOS MR, however, our results 
remained ambiguous even after long integration times. Thus, 
we moved to a more detailed procedure for deciding which 
mass shells were ejecta and which composed the remnant 
(see Fig.[3J. We stored the location and local energy density 
of each grid zone throughout the run. We then plotted the 
location of each mass element as a function of time, using 
the sign of the local energy density to color code the lines. 
During the atmospheric motions some layers do work on 
other layers; the color coding shows changes from bound to 
unbound (and vice- versa). These plots proved to be helpful, 
illuminating the transient identities of bound atmospheres, 
marginally bound gas, and low energy ejecta. We adopted 
the following criterion for ending the calculation: when all 
outer shells had positive energy density and the number of 
intermediate shells with local energy density still changing 
sign was small - less than a couple of percent of the total 
mass. An example is shown in Fig. [3] where the apparent 
bifurcation point between bound and unbound material is 
marked on the far right. Note the diminishing amount of 
mass ejected with each stellar oscillation. We are confident 
of this prediction because it was borne out in all cases where 
the code was run much longer, and, hence, closer to the point 
of the remnant's return to hydrodynamic stability. 

These plots illustrate two distinct ways in which shells 
are ejected in successful explosions that are not yet large 
enough to destroy completely the stellar model: 1) Some 
shells are lost in the initial shock wave; these gain tremen- 
dous kinetic energies. The amount of mass lost in this way is 
typically < 12 per cent. 2) The rest of the ejecta are expelled 
by the ringing down of the remnant. 



Figure 3. A typical graph of the motion of grid zones in time for 
an exploding stellar model, in this case an rh = 100 polytrope of 
index n = 3 with an explosion energy of 90 per cent of the stellar 
model's binding energy. Each line tracks a representative mass 
shell. The Lagrangian mass intervals vary: lines in the ejected 
region and outer parts of the remnant — which represent escap- 
ing mass (large radii) and the uppermost parts of the cooling, 
bouncing remnant (the blue lines) - are at intervals of 0.5 — 1 per 
cent of the total mass. In the inner part of the remnant, each line 
represents approximately 10 per cent of the total mass. Where 
lines are red, the local energy density is positive; where blue, 
negative. The bifurcation point separating the remnant from the 
ejecta is marked. Radial distances are given in units of the initial 
stellar radius; times are given in dimensionless units defined by 
'/(RL r /GM Jt „ r ) 1/2 - 

We next investigated the extreme limits: total disrup- 
tion explosions and failed explosions (no ejected mass). Total 
disruptions were relatively easy to recognize: all grid zones 
acquired positive energy in the first pass of the shock wave 
from the explosion. For explosion energies near the thresh- 
old for total stellar disruption, however, we found that there 
was typically a range of energies where a large portion of 
the stellar envelope was ejected, while leaving an extended 
remnant undergoing very slow, long-scale oscillations. The 
precise mass loss in most of these cases was impossible to 
determine, though it was always ~ 50 per cent. Artificial 
viscosity eventually damps these oscillations, but it takes a 
long time to do so. Total disruption happens abruptly, with 
every stellar model studied going from the ~ 50 per cent 
mass loss oscillatory state to 100 per cent mass loss with 
only a small increase in explosion energy. We were able to 
pin down the width of the transition from remnant to total 
disruption as function of explosion energy to ~ 5 — 10 per 
cent in the stellar model's binding energy. 

The failed explosion regime was computationally eas- 
ier to study. Failed explosions produced no unbound shells. 
The results can be understood in terms of the varying shock 
speed. Strong shocks slowed as they plowed through the 
dense core of the stellar model, then accelerated when they 
reached the steep density gradient of the outer regions of 
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Figure 4. A plot of the v gas / \J v^ sc + ° 2 sn d' * n ^his case f° r a 
stellar model with n = 1.5 , fa = 0.55, and an explosion energy 
equal to 15 per cent of the stellar model's binding energy. The 
choppiness in the plot is due to non-physical effects in the deter- 
mination of the exact location of the shock. Note the deceleration 
through the bulk of the stellar model, with only the very outer- 
most shells reaching escape velocity as the shock accelerates in 
the falling density profile near the edge. 

the polytropes. In failed explosions the shock velocity fell 
below the sound speed in the middle region and/or failed to 
accelerate up to the local escape speed in the outer region. 
A plot of the process is included in Fig^] In the figure, we 
plot v shock / yj vl sc + c 2 and - where v e3C is the escape velocity 
for the initial stellar model and c sn d is the local sound speed 
- illustrating the falling shock Mach number in the core and 
the reacceleration to velocities allowing escape by the large, 
negative, outward-going density gradient. The physics of this 
variable shock velocity, including spheri cal curvature effects 
and photon loss, are discussed in detail in M atzner fc McKed 
(1999). 

4.4 Explosions without radiation pressure 

For the first round of explosions, we used EOS M (no radi- 
ation pressure). For these calculations, the solution's inde- 
pendent parameters are -Ebiast / -Ebmd and the stellar model's 
polytropic index n. 

We have included two figures summarizing the mass loss 
results. In the first, Fig- El the explosions are compared with 
each other, showing great similarity among the models. In 
Fig. |S] we have separated each polytrope into its own win- 
dow to c ompare its mass loss cur ve to a couple of "rules of 
thumb" jFrver fc KaloeeralEoOll) . The first line, the dash- 
dotted curve, is the simplest such rule. It represents the mass 
loss if 100 per cent of the explosion energy were distributed 
in such a way as to eject as many of the outer shells as 
possible, while leaving untouched those parts of the stellar 
model which remain bound. This is, of course, physically 
impossible, but it does provide an upper bound on mass 



Figure 5. This figure summarizes the mass loss percentages re- 
sulting from explosions in polytropes of four different indices with- 
out radiation. 

loss. The dashed curve r epresents the actual choice made by 
iFrver fc Kaloeeral <l200ll) . which essentially splits the explo- 
sion energy budget in two, giving 50 per cent to unbinding 
the stellar model directly, and 50 per cent to heating the 
remnant and to accelerating the ejecta. This version gives 
results that are much closer to our numerical calculation, 
but overestimates mass loss in low energy explosions and 
also overestimates the amount of energy required to unbind 
the stellar model completely. 

4.5 Explosions with radiation pressure 

For the second round of explosions, we used the hydrody- 
namics code with EOS MR (matter and radiation pressure). 
The parameter space now included three variables: explo- 
sion energy, polytropic index, and M s t ar - We chose three 
stellar masses, rh = 10, 100, and 1000, and then surveyed 
the range of polytropic indices for convectively stable mod- 
els. For each stellar model, we varied explosion energy from 
cases of failed explosions to total stellar disruption. The re- 
sults of these explosions are summarized in Figs. |7| |H| and 
|U] Because of the difficulty of precisely determining the line 
of bifurcation between remnant and ejecta even in plots like 
Fig. [3] error bars have been included. They represent the 
range within the stellar model where the bifurcation point 
may occur. This range of uncertainty was determined by 
finding the region of the stellar model containing either out- 
going zones with negative local energy density, located at 
several stellar radii, or infalling zones with positive local en- 
ergy density. We considered these to have ambiguous fates. 
The number of these zones is always small, comprising at 
the most a couple of percent of the stellar mass. Bold dots 
are placed at the midpoint of this range. 

The mass-loss curves are remarkably uniform, especially 
given the wide variation in binding energy among the stellar 
models that we studied (more details on binding energy are 
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Figure 6. Mass loss curves are compared with two simple as- 
sumptions relating explosion energy and mass loss. The dashed 
curve represents the most efficient possible application of the ex- 
plosion energy to mass loss. The dash-dotted curve represents the 
more physically reasonable assumption that 50 per cent of the en- 
ergy goes into unbinding part of the stellar model, and 50 per cent 
goes into both heating the remnant and to net kinetic energy for 
the ejecta. 



Figure 8. This figure summarizes the mass loss percentages for 
explosions in fa = 100 stellar models. Error bars represent the un- 
certainty in the determination of the remnant - ejecta bifurcation 
point. 
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Figure 7. This figure summarizes the mass loss percentages for 
explosions in fa = 10 stellar models. Error bars represent the un- 
certainty in the determination of the remnant - ejecta bifurcation 
point. 
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Figure 9. This figure summarizes the mass loss percentages for 
explosions in fa = 1000 stellar models. Error bars represent the 
uncertainty in the determination of the remnant - ejecta bifurca- 
tion point. 



given in Appendix B). The uniformity in the shapes of the 
mass-loss curves found allows them to be described accu- 
rately by a fitting formula. The form that best fits the data 



M f e < e D 

100 x— 12!1 = J A(e blast - e D ) 2 e D < e < e f (7) 

aw y 100 ej<e 

where euast (e Q , e/) is blast energy (minimum blast energy 
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Figure 10. A figure comparing the results of the fitting formula, 
Eqn. \7\ to the data, in this case for an n = 3, rh = 100 stellar 
model. 



Table 2. Summary of fitting parameters for various explosion 
scenarios. e a = 100 X -Eo/£bind> e f = 100 X -E//-Ebind 



rh A(xl0 -3 ) 



e f 



2.0 


10 


3.57 


10 


110 


2.5 




3.37 


11 


120 


3.0 




3.47 


10 


120 


4.0 




3.67 


10 


130 


2.5 


100 


2.89 


5 


120 


3.0 




3.50 


10 


130 


4.0 




2.95 


11 


140 


3.0 


1000 


3.74 


20 


110 


4.0 




2.44 


10 


130 



to cause mass loss, maximum blast energy to leave a bound 
core) measured as a percent of binding energy (i.e. ebiast = 
100 x -Ebiast /Sbmd, etc.) and A is the fitting parameter. A 
sample comparison between the fits and the numerical data 
is shown graphically in Fig. 1101 the parameters describing 
each model's fit are contained in Table [5] 

Notice that in all cases, total disruption occurs near 
or slightly above the original stellar binding energy, i.e. at 
Ebiast = (e//100)Ebind ^ Ebind- The transition to total dis- 
ruption is also very abrupt. Just below e/, the amount of 
mass lost is below 50 per cent in all cases. 



None of our models included a compact object at the cen- 
tre and one naturally wonders whether the later mass loss 
might be sensitive to the treatment of the centre of the stel- 
lar model. 

To test the effect of the inner regions on the mass loss, 
we considered two different ways of modifying the inner 
boundary condition. In one, we placed a hard, reflecting 
sphere at a fixed, very small spatial radius; in the other, we 
fixed a perfectly absorbing boundary at a given radius. We 
ran several cases: polytropes of indices n between 1.5 and 4, 
with a variety of values of rh. For the reflecting sphere, the 
size was set to that of the innermost zone; for the absorb- 
ing sphere, the absorbing boundary varied between 2 and 30 
percent of the mass in Lagrangian coordinates. 

The mass loss computed with these altered inner bound- 
ary conditions differed little from what was found in our 
earlier survey. The changes to the mass loss were unnotice- 
able even when quite a large inner region - as much as 15 
- 20 per cent, in mass - of the stellar model was allowed 
to become pressure-free and perfectly absorbing. This in- 
sensitivity is linked to the fact that at least 50 per cent of 
the star remains bound if the star survives the explosion. In 
those cases, on the one hand, the trading of energy among 
the outer shells determines the amount of mass lost; the in- 
ner half of the stellar model, whose dynamics are sensitive 
to the treatment of the inner zones, is never in danger of be- 
ing lost. In stellar models that are totally dispersed, on the 
other hand, the explosion is so violent that differences in the 
inner boundary conditions have little effect, since no portion 
of the stellar model ever falls back to experience them. 

To gain a better qualitative understanding of why the 
mass loss results are so insensitive to the inner boundary 
condition, we examined more closely the slowly-varying dis- 
turbances propagating in our remnant. The energy density 
in a small- wavelength oscillatory mode is proportional to lo- 
cal density times the square of local gas displacement (e.g. 
IChristensen-Dalsgaardl 1 20 03) . We computed this combina- 
tion - as well as local velocity squared times local density - 
as a function of radius at various times in our models, using 
each Lagrangian zone's displacement from its final, at-rest 
location. We found that both quantities always go to zero 
at the remnant core, implying that there is little energy flux 
into or out of this region. Therefore, information about the 
altered inner boundary conditions is not effectively transmit- 
ted to the oscillating outer mass shells. The inner boundary 
condition plays no obvious role in the late-time mass loss so 
long as there is a sufficient amount of buffering gas between 
the inner core and the outside of the star. In our numeri- 
cal tests, we found that it was necessary to alter almost the 
entire inner quarter of the stellar model's mass before the 
mass loss was significantly affected. 



5 SENSITIVITY TO INNER BOUNDARY 
CONDITIONS 

Our previous results show that mass loss occurs in two dis- 
tinct phases when the stellar model is not completely un- 
bound. The initial shock wave drives off some mass immedi- 
ately while the ringing down of the remnant expels loosely 
bound outer shells over a period of several dynamical times. 



6 CONCLUSION 

Here, we have confined ourselves to a rather idealized prob- 
lem, partially disruptive explosions of simplified stellar mod- 
els whose density profiles are solutions to the Lane-Emden 
equation, but with varying ratios of radiation to matter 
pressure. In this way, we have been able to survey models 
in a well-defined parameter space fairly extensively, trading 
the concreteness of real stellar models for the flexibility of 
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parametrized ideal models. Specifically, we have determined 
the fraction of the original stellar mass ejected as a func- 
tion of explosion energy in polytropes of index n = 1.5, 
2, 3, and 4 in calculations without radiation pressure; we 
also explored the mass loss fractions for stellar models of 
in — 10, 100, and 1000 over a range of n from 2.0 to 4.0. 
Our results suggest that the mass loss is remarkably uni- 
form, even among widely varying stellar density profiles and 
over an enormous range in stellar masses. We have provided 
a simple, parametrized formula for the fractional mass loss 
as a function of explosion energy for a range of values of n 
and the ratio of radiation to matter pressures; see Eq. JJJ 
and Table H 

One striking feature of all the models we tested was that 
the mass loss fraction as a function of explosion energy ap- 
pears to be discontinuous at around 50 per cent mass loss, 
with a small (few per cent) difference in explosion energy 
separating stellar models which lose half their mass from to- 
tally disrupted stellar models. Because our simulations did 
not include formation of a compact remnant at the centre, 
this result cannot be taken as a concrete demonstration that 
the observation of a black hole of mass M demands a pro- 
genitor whose mass was less than about twice as large. 

The abrupt transition between moderate (i.e. < 50 per 
cent) and total disruption found here for wide classes of ini- 
tial models is also seen i n modelling of explosions in sets of 
specific progenitors (e. g. IWooslev fc Weaverlll995l Table 3; 
iMacFadven et al.ll200ll Table 1). Thus, we conjecture that 
even when a compact central remnant is included, the results 
divide into two separate cases depending on whether the ex- 
plosion energy is above or below, approximately, the critical 
value found here for complete disruption. It is significant 
that this bifurcation effect occurs, as it does in our simula- 
tion, at the most basic, hydrodynamic level, implying that 
its appearance in more sophisticated models rests chiefly 
on these underlying physics. For explosion energies below 
this critical value, there is a sharp transition between mod- 
est (i.e. < 50 per cent) mass loss and total disruption apart 
from the compact remnant. For explosion energies above the 
cutoff, either a black hole or neutron star may form. How- 
ever, in this case, we expect much smaller fallback masses, 
generally only a few tenths of Mq or less, primarily caused 
by reverse shock propagation through the core, and the con- 
sequ ent deceleration of a small am ount of outgoing matter 
fe.p. lWooslevlll988l: IChevalierlll989T) . 

We should emphasize that our conjectures here may 
be significantly modified by the inclusion of more realis- 
tic physics. Spherically-symmetric polytropic stellar models 
are not real stars; neither are all the complexities of stel- 
lar hydrodynamics captured by simple equations of state. 
One important component of realistic stars that our models 
lack is the sequence of density discontinuities that emerge 
in evolved stars, which we would expect to cause reflected 
reverse shocks that could change the mass ejection dynam- 
ics considerably. Furthermore, in real supernovae there are 
other sources of energy besides the initial explosion - such as 
decaying radioactive nuclei, central pulsars or jets - which 
continue to add energy to the supernova system at late 
times, possibly giving the boost needed to eject lightly- 
bound gas at the surface of the remnant. Our imposition 
of spherical symmetry is also unrealistic: real explosion are 
likely to be anisotropic, and stellar rotation, convection, 



and magnetic fields are expected to have effects that a one- 
dimensional calculation cannot hope to model. Despite these 
caveats, we believe our chief results to be robust. 

Supernovae frequently occur in binary systems. In such 
systems, when there are very energetic explosions and little 
mass fallback, we expect little mass contamination of the 
atmosphere of the binary companion. The outgoing regions 
of the progenitor intercepted by the companion are not cap- 
tured; indeed the outer layers of the companion are stripped 
and ablated by the ejecta. On the other hand, when the 
explosion is weak and substantial mass fallback occurs, pro- 
genitor material may fall back onto the companion, polluting 
its atmosphere. In the latter cases, we would then infer that 
a remnant of mass M was most likely derived from a pro- 
genitor with mass less than ~ 2M. Thus, in systems like 
Nova Scorpii that show evidence for bl ack hole formation in 
a supernova (e.g. Ilsraelian et al.lll999l) . we conjecture that 
mass of the pre-explosion star was, in fact, less than twice 
the present mass inferred for the black hole remnant (which 
also has accreted matter since forming, presumably). This 
may have implicatio ns for the dynamics of such systems (e.g. 
iMirabel et alJ 120021 1. We caution, though, that our results 
may be altered somewhat in more refined models. Further 
studies are underway to include a compact central remnant, 
density jumps (expected as a consequence of compositional 
inhomogeneity) , rotation and explosion asymmetries. These 
new calculations will continue, in the same spirit as those re- 
ported here, to employ the simplest explosion models needed 
to reveal the underlying physical consequences of the various 
refinements, and to allow a survey of the hydrodynamics of 
a large range of explosion models. 
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APPENDIX A: DIFFERENCE EQUATIONS 

The Lax-Wendroff difference equations for the equations of 
hydrodynamics in one dimension with spherical symmetry 
are as follows. Note that the pressure in the equation for 
advancing energy must be solved for using the equation of 
state to make this set of difference equations explicit rather 
than implicit. In the difference equations, n represents time 
steps, while j represents spatial steps. The equations are non- 
dimensionalized simply, with each variable scaled to order 
unity for the initial conditions in all calculations we have 
done. The one remaining constant, p , with units of density, 
sets the overall scale of the system studied. The variable 
R records the position of each shell. Comparing each shell's 
current position, R, with r, a static, reference grid, allows the 
gas's local density to be calculated. The remaining variables 
are interdependent. The equation for moving grid zones is: 



At 1 
The conservation of momentum equation is: 



At 



Po 



Ar 



The conservation of mass equation is: 



n+l _ 
Pj + i/2 — P° 



(7?^+ 1 1 )3-(^+ 1 )3- 



(Al) 



(A2) 



(A3) 



The First Law of Thermodynamics is: 



U 



n + l 

j + 1/2 



3 + 1/2 



-I- n" 
/2 ^ Pj + l/2 



n+l n n 
P 3 + 1/2 P 3 + 1/2 



(A4) 



Where U = internal energy / mass. The acceleration of the 
innermost shell is determined by treating its volume as filled 
with a gas of uniform pressure so that the shell's equation 
of motion is: 

dv 



Thinner — ^7f(,Pinner Pouter) 



n+l 
U 



nner Pouter)- 



(A5) 



The pressure within the inner sphere varies adiabatically as 
the shell moves, i.e., 



Pinner(t) — Po 



Vo 

V(t) 



(A6) 



These equations are completed by some equation of 
state, 

n+l 



n+l _ ftrjn + l- „n+l \ 
Pj + 1/2 - J( U j + l/2'Pj+l/2>- 



(A7) 



If this equation of state can be algebraically solved, the full 
set of equations is explicit; if it cannot be solved, then an im- 
plicit step and numerical root-finding procedure is required 
to advance the grid. The advancement of the grid proceeds 
as follows. 1) Using the conservation of momentum, the new 
gas velocities are set throughout the system. 2) Boundary 
conditions are applied. 3) The shell position, R, is advanced 
according to the new gas velocities. 4) R is then used to set 
the density throughout the system. 5) Two possibilities: if 
the equation of state is explicitly soluble, then the internal 
energy of the gas is determined. If not, then the pressure 
and energy equations must be stated in terms of the tem- 
perature and then solved, together with the First Law, nu- 
merically - three equations for three variables, p, U, and T. 
The above prescription must be modified slightly to accom- 
modate shock fitting. To this end, we introduce an artificial 
viscous pressure, q, given by the differenced form, 



f 2a 2 [(8u)™ +1/2 ] 2 
<h+l/2 = \ 1 l»l+l,2 + 1 lf>l+t / 2 





if (Su) 
if (Su) 



3 + 1/2 



< 

> 



(A8) 



) 3 + l/2 

Note the parameter, a, which controls how widely the vis- 
cous pressure spreads the shock. Optimal values are 1.5 < 
a < 2.0, which spread the shock over 3-10 zones. This artifi- 
cial viscous pressure is added to the regular gas pressure in 
the above equations as follows: In the conservation of mo- 
mentum equation, 



( s p)7 -> ( s p)] + ( s i)7 

and in the energy conservation equation, 

j+1/2 P j + 1/2 ^ Pj + 1/2 



„n+l , • 
Pj + 1/2 + P 



+ a n+1 

+ 13 + 1/2- 



(A9) 



(MO) 



2 2 
When advancing the grid with artificial viscous pressure, 
the artificial viscosity term, q, is advanced before the energy 
equation, step 5 in the previous description. 



APPENDIX B: SETTING UP MASS SHELLS IN 
A POLYTROPE 

In the initial configuration whose mass and radius are M* 
and R*, define a mass coordinate rh — M/M* so that the 
shell is at radius r(m) = R(M)/R+. 

Let the pressure and density be P(M) = 
P{rh)(GM?/Ri) and p(M) = p(m)M*/R'l, respectively. At 
the centre of the stellar model, we can find P(0) = fp(n) 
and p(0) = f P (n), where n is the polytropic index. At any 
other point in the stellar model, we have p(M) = p(0)8 n 
and P(M) = P(0)8 n+1 . Thus, we have 



P(m) = f P (n)[6(m)Y 



p(m) = f P (n)[8(m)Y 



(Bl) 



Usually, we specify the Lane-Emden function as a function 
of a dimensionlcss radius. Getting r(m) then requires a little 
bit of work. The method is this: the mass inside a physical 
radius R is 



7" 

^0 



M(R) = 477 / dr r 1 p(r) = 47rr s d calc / dx x 2 [8(x)] n , (B2) 
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where r sca i c is the radius scale. Divide by the total mass to 
get the conversion equation 



xo(n)r(rh) 



dx x 2 [8{x)] n 



r; o(n) dx x 2 lew /; ow dx x 2 i8(x)\ 



x (n) 



(B3) 



where 8(xo(n)) — and we have used the fact that the scaled 
Lane-Emden radius variable x = xo(n)r(m). This equation 
can be inverted numerically to get f(m), and we can use the 
result to evaluate the pressure and density: 

P(rh) = f P (n)[e(x (n)f(m))r +1 
P(m) = f P (n){8(x (n)r(m))] n . 

We are interested in polytropic models where the pres- 
sure is supplied by a mixture of radiation and a nonrelativis- 
tic gas. The total pressure in physical units is then 

P = — h -^ a T A = -Pgas + -Pradiation , (B4) 

jji O 

where p is the mass per nonrelativistic particle in the stellar 
model. We define s 



aT A p 



(B5) 



-Pgas 3pk 

which varies throughout the stellar model. At any point in 
the stellar model, we can use this to eliminate temperature 
in favor of s: 

T= (M 1/3 ^ 

V ay, / 

p 



pkT(l+s) 



From this it follows that 

l+l/n 
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P(0) 
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(p(o)) 



4/3 
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1/3 



1 + S 



[8((xo(n)r(rh))] in/3 



1 + s(0) 

1/3 



1 + s 



We can use this to solve for s(m) via 

[s(m)] 1/:i [l+s(m)] = [s(0)] 1/3 [l+ S (0)] [8 (lot^fW)] 1 ""/ 3 . 

The thermal or internal energy density inside the stellar 
model is, in physical units, 



pU 



2p 



2p 



- !(?r©" s <-> 



compare this with the pressure 
Define 7(rh) by P = (7(m) — l)pU, to find that 



P = 



(t) (1 + s) - 



pU 
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let U = (GM*/ Ri,)U to find that 



2/ p (n) 



1 + 2s(m) 



1 + s(m) 



(B9) 



This completes the setup of the initial conditions of the poly- 
trope. 

We will also need to have the total energy of the stel- 
lar model because we want to choose the blast energy as a 
fraction of the binding energy. The gravitational energy of 
a polytrope of index n is 

3GM 2 

-C^grav 



(5-n)i?* ■ 
The internal energy of the stellar model is 
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1 + s(m) 



This is enough to get the total energy, but we can be a slight 
bit more elegant by using the virial theorem, 



-Bgrav = 3 I dM 



3f P {n)GMl f 1 
= f P {n)R+ J d(x (n)r(m)) , (Bll) 



from which we find that 
3GM 2 



^tot 
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J dm 8(x (n)f(m)) 

If we define E tot = -E tot [3GM 2 /2(5 - n)R+] we see that 
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APPENDIX C: DYNAMICAL EQUATIONS 

The equation of motion for a mass shell is 



d 2 R(M,t) 
dt 2 



= -4TvR 2 (M,t) 
GM 



dP(M,t) 
dM 



+ a viBC (M,t) , 



(CI) 



R 2 {M,t) 

where a v j sc (M, t) is the viscous acceleration (which we in- 
clude using a prescribed artificial viscosity) . Introducing our 
nondimensional radius, pressure and mass implies 

d 2 f(m,t) 
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GM* 
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-47rf (m, t) 
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dm 
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f 2 (m, t) 

Define a dimensionless time by t = (i? 3 /GM*) 1 ' 2 r; then 
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! (m, t) 



+ a viBC (m, t) 
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where the viscous acceleration is defined by a V isc(M, t) = 
a v i BC (m,T)(GM*/Rl). Since we shall actually want equa- 
tions that are first order in time, we note that the radial 
velocity is 



dR(M, t) 
dt 



and therefore 

df(rh, t) 

cV 
dv(m, t) 



/ GA^ \ 1/2 df(m,r) 
\~~Rr) ~d~? 

/GM,\ 1/2 



v(m, t) 



(C4) 



v(m, r) 



dr 



, ~2, - \ dP(m, t) 
-vkt (m, r) 



2 (m,r) 



dm 
+ a visc (m,T) 



(C5) 



From the first law of thermodynamics, we get that 
dU(M,t) . %r . n ,, f .S [ f 

where q v i sc (M, t) is the viscous heating, which we take to be 
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g visc (M,t) = 
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where a is a constant with units of length. Introducing the 
same nondimensional variables as in the polytrope setup we 
find that 

dU(m,r) ( 2 „. „ , (dv\ 2 A . „ 

= - ap(m,T)[- J +P(m,T) ] > 
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where the viscous heating is defined by g v j sc 
9viac[(GM Vr ) 3/ ' 2 /i?^ ]. We can replace the pressure by 



\l + 2sj 
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to rewrite the first law in the form 
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We can close the loop using (7 oc p 1//3 s 1/ ^ 3 (l + 2s) to get 



(C9) 



U(m,T) f p(m,T)s(m,T) 



U(m,0) \p(m,0)s(m,0) 



1/3 



1 + 2s(ra,r) 



1 + 2s(m,0) 



(CIO) 



We could use this to find an equation for s(m, t) explicitly; 
then 



d_ 

Or 



s(rh, t)< 



,8s(m,r) 



p(m,r) 



3g v isc(fn,r)[l + 2s(m,r)} 
U(m,r) 



p{m,T) 



(Cll) 



with U (m, t) evaluated using Eq. IClOt . Finally, the equa- 
tion of mass conservation can be written as 



df(m, t) 



p(M,t) 
1 

p(fh,T) 



4:Tvf (m, t) 



dm 



(C12) 



If we wish, we can define a new variable V(m,r) 
1/ p(m,r), and rewrite the last few equations as 



V{m,r) 

U{m,r) 
U(m,Q) 
dU(m, t) 



— 47rf 2 (m, r) 



df(m, t) 



dm 

V(m, 0)s(m, t) 



1/3 



V(m, r)s(m, 0) 



fdvy 



1 + 2s(m,r) 



1 + 2s(m,0) 



+ 



V(m, t) 
2U(m,T)[l + s(m,r)] 



dV(m,T) 



(C13) 



3V(m,r)[l + 2s(m,r)]J dr 

Eqs. dC5l . and either Eqs. EU, dcTol and iC12t or Eqs. 
l|C13|l . with the initial conditions set up in the previous sec- 
tions, can now be cast into finite difference form, with a 
suitable specification of the artificial viscous force and heat- 
ing. 

For all calculations done after the code was tested, we 
also included a Newtonian Gravitation force per unit mass 
via 

GM{R) 



Fgrav — 2 , 

or in difference form, 
Gpof (r?) 3 



F n — _ 
1 j 



(C14) 



(C15) 



This force was added to the conservation of momentum 
equation, the equation used to set gas velocities. Finally, 
the code self-checks by calculating total energy and momen- 
tum to ensure that these are conserved. For energy, the sum 
of the local energy in each zone is calculated first via 



Skin 



+ E ti 



E 

i^_zones 



(^u 2 + U t )AMi. 



(C16) 



The gravitational potential energy is then calculated via 

nclosed 



i^zones 



Ri 



-AM;. 



(C17) 



and the two energies are added and recorded as the current 
total energy in the system. Conservation of momentum is 
also checked though a simple summation: 

Ptot = u zone AMi. (C18) 

zones 

Finally, the algebraic equation used to determine the total 
local energy per unit mass of each zone - the quantity used 
to determine if a zone was bound or unbound - was 

potential 



E 



■prbherm | ^kinetic. _j_ jrrp 



AM, 



2 U ? + V i 



where ARi = R t — Ri-i. 



AM j 

i=l 



F? rav ARi 



(C19) 
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APPENDIX D: SEDOV SOLUTION 

For the analytic solution to the Sedov problem, we rederived 
the solution given in Landau and Lifshitz's Fluid Mechanics, 
thereby finding the correction to an error (in an exponent) 
that is mentioned in Shu's Gas Dynamics. Rather than re- 
reporting the entire result, we simply note that the correc- 
tion comes in the exponents of the equation for p: 

-((^) (-§"))*' 

where 

2 



